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Abstract — The finite-difference time-domain (FDTD) method 
is applied for modelling of wire media as artificial dielectrics. 
Both frequency dispersion and spatial dispersion effects in wire 
media are taken into account using the auxiliary differential 
equation (ADE) method. According to the authors' knowledge, 
this is the first time when the spatial dispersion effect is consid- 
ered in the FDTD modelling. The stability of developed spatially 
dispersive FDTD formulations is analysed through the use of von 
Neumann method combined with the Routh-Hurwitz criterion. 
The results show that the conventional stability Courant limit is 
preserved using standard discretisation scheme for wire media 
modelling. Flat sub-wavelength lenses formed by wire media 
are chosen for validation of proposed spatially dispersive FDTD 
formulation. Results of the simulations demonstrate excellent sub- 
wavelength imaging capability of the wire medium slabs. The 
size of the simulation domain is significantly reduced using the 
modified perfectly matched layer (MPML) which can be placed 
in close vicinity of the wire medium. It is demonstrated that 
the reflections from the MPML-wire medium interface are less 
than -70 dB, that lead to dramatic improvement of convergence 
compared to conventional simulations. 



I. Introduction 

The wire medium is an artificial material formed by a 
regular lattice of ideally conducting wires (see Fig. Q. The 
radii of wires are assumed to be small compared to the lattice 
periods and the wavelength. The wire medium has been known 
for a long time [l]-[3] as an artificial dielectric with plasma- 
like frequency dependent permittivity, but only recently it was 
shown that this dielectric is non-local and possesses also strong 
spatial dispersion even at very low frequencies [4]. Following 
[4] the wire medium can be described (if lattice periods are 
much smaller than the wavelength) as a uniaxial dielectric with 
both frequency and spatially dependent effective permittivity: 



e = e{k, ga;)xx+yy+zz, e{k, q.^) = £o 1 






fc2 



(1) 



where k^ = wg/c is the wave number corresponding to the 
plasma frequency ojq, k = lu/c is the wave number of free 
space, c is the speed of light, and q^ is the component of wave 
vector q along the wires. The dependence of permittivity Q 
on q^ represents the spatial dispersion effect which was not 
taken into account in the conventional local uniaxial model of 
wire medium [l]-[3]. 

The plasma frequency of wire medium depends on the 
lattice periods a and b, and on the radius of wires r [5]: 
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Fig. 1. The geometry of the wire medium: a rectangular lattice of parallel 
ideally conducting thin wires 
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For the commonly used case of the square grid (a = b), 
F(l) =0.5275. 

The finite-difference time-domain (FDTD) method has been 
widely used for modelling of transient wave propagation in 
frequency dispersive and non-dispersive media [6]. The exist- 
ing frequency dispersive FDTD methods can be categorised 
into three types: the recursive convolution (RC) method [7], 
the auxiliary differential equation (ADE) method [8] and 
the Z-transform method [9]. The first frequency dispersive 
FDTD formulation was developed by Luebbers et al. for 
the modelling of Debye media [7] using a RC scheme by 
relating the electric flux density to the electric field through 
a convolution integral, and then discretising the integral as 
a running sum. The RC approach was also extended for the 
study of wave propagation in a Drude material [10], M-th order 
dispersive media [11], an anisotropic magneto-active plasma 
[12], fenTte material [13], and the bi-isotropic/chiral media 
[14]-[16]. The ADE method was first used by Kashiwa and co- 
workers [17], [18] in 1990 for modelling of Debye and Lorentz 
media. Taflove [19] soon extended this model to include effects 
for nonlinear dispersive media. Independently, Gandhi et al. 
proposed the ADE method for treating M-th order dispersive 
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media [20], [21]. The other dispersive formulation based on Z 
transforms was proposed by Sullivan [9] and then extended to 
treat nonlinear optical phenomena [22] and chiral media [23]. 
In [24], Feise et al. compared the ADE and Z transform meth- 
ods and applied the pseudo-spectral time-domain technique for 
the modelling of backward-wave metamaterials to avoid the 
numerical artifact due to the staggered grid in FDTD. Lu et 
al. used the effective medium dispersive FDTD method for 
modelling of the layered left-handed metamaterial [25] and 
demonstrated the zero phase-delay transmission phenomenon. 
Recently, Lee et al. applied the piecewise linear RC (PLRC) 
method through an effective medium approach for modelling 
of the left-handed materials using the similarity of its effective 
permittivity and permeability functions to Lorentz material 
model [26]. 

The simulation of wire medium can be performed either by 
modelling the actual structures i.e. parallel wires, or through 
the effective medium approach. However, in order to accu- 
rately model thin wires in FDTD, either enough fine mesh 
is required or the conformal FDTD [27] needs to be used. In 
this paper, the effective medium approach is used and the wire 
medium is modelled as a dielectric material with permittivity 
of the form Q. Due to the similarity of the frequency and 
spatial dispersion effects in ([ij with Drude material model, 
the ADE method can be directly applied for the spatially 
dispersive FDTD modelling. 

II. Dispersive FDTD Formulations 

Using the Eq. Q the wire media can be modelled in FDTD 
as a frequency and spatially dispersive dielectric. In order 
to take into account the dispersive properties of materials 
in FDTD modelling, the electric flux density is introduced 
into the standard FDTD updating equations. Since the x- 
component of electric flux density Dx{u!, Qx) is related to the 
cc-component of the electric field intensity Ex{iLJ,qx) in the 
spectral (frequency-wave vector) domain as 



Dx{uj) = e{uj,qx)Ex{LL>), 



one can write that 



(e - ql) D, + {ql -e + kl) eoE^ = 0. 



(4) 



(5) 



This equation allows to obtain the constitutive relation in the 
time-space domain in the following form: 



\dx^ 



D,, 
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dx^ 



k' eoE, = 0, 



(6) 
using inverse Fourier transformation and the following rules: 

The Eq. (|6} relates only a;-components of the electric flux 
density and field intensity. The permittivity in both y- and 
z-directions is the same as in free space since the wires are 
assumed to be thin. 

The FDTD simulation domain is represented by an equally 
spaced three-dimensional (3-D) grid with periods A^, Ay 
and A^ along x-, y- and z-directions, respectively. The time 



step is Aj. For discretisation of (|5), we use the central finite 
difference operators in time (6^) and space (5^) as well as the 
central average operator with respect to time (/ij): 
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where the operators 51 , 5'^ and ixj are defined as in [28]: 
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Here F represents the field components; mx,my,mz are 
indices corresponding to a certain discretisation point in the 
FDTD domain, and n is the number of the time steps. The 
discretised Eq. (|6) reads 
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Note that in (|8j, the discretisation of the term fcg of (|6j 
is performed using the central average operator ^^ in order 
to guarantee the improved stability. The stability of different 
discretisation schemes is analysed in details in the next section. 
The Eq. (|8} can be written as 
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Therefore, the updating equation for E, in terms of E, and 
Dx at previous time steps is as follows: 
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with the coefficients given by 
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Fig. 2. The layout of the computation domain for two-dimensional FDTD 
simulations 
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In the spatially dispersive FDTD modelling of wire medium, 
the calculations of D from the magnetic field intensity H, and 
H from E are performed using Yee's standard FDTD equations 
[6], while Ex is calculated from Dx using \\Q\ and Ey — 
e^^Dy, Ez ~ Eq^Dz- Note that in ( I10> . the central difference 
approximations in time (for frequency dispersion) for both Dx 
and Ex are used at position {mx,my,mz), and the central 
difference approximations in space (for spatial dispersion) are 
used at the time step n in order to update Ex at time step 
n + 1. Therefore the storage of Dx and Ex at two previous 
time steps are required. 

At the free space-wire medium interfaces along x-direction, 
the updating equation ( I10> includes Dx and Ex at previous 
time step in both free space and wire medium. Outside the 
region of wire medium, the updating equation (^J reduces to 
the equation relating Dx and E^ in free space. 

The spatially dispersive FDTD method has been imple- 
mented in a two-dimensional case and used for modelling 
of sub-wavelength imaging provided by a finite sized slab of 
wire medium [29]. The simulated geometry is illustrated in 
Fig. 12] Different types of sources are chosen in simulations 
including a magnetic point source and three equally spaced 
magnetic point sources in order to demonstrate the sub- 
wavelength imaging capability of the device. The simulation 
results are provided in section V, while the stability and 
numerical dispersion relation for general 3-D case are analysed 
in the following section. 

III. Stability and Numerical Dispersion Analysis 

A. Stability 

Previous stability analysis of dispersive FDTD schemes is 
performed using the von Neumann method and numerical root 



searching [30]. In this paper, the stability of the proposed spa- 
tially dispersive FDTD method is analysed using the method 
combining the von Neumann method with the Routh-Hurwitz 
criterion as introduced in [31]. The von Neumann method 
establishes that, for a finite-difference scheme to be stable, all 
the roots Zi of the stability polynomial S{Z) must be inside 
of the unit circle in the Z-plane (i.e. \Zi\ < 1 V i)> where the 
complex variable Z corresponds to the growth factor of the 
error and is often called the amplification factor [31]. 

The wire medium is an uniaxial material where the di- 
vergence of electric field inside wire medium is non-zero 
(V • E 7^ 0). Therefore, to analyse numerical stability of the 
proposed spatially dispersive FDTD method, we must start 
directly with the Maxwell's equations instead of the wave 
equation as was done in [31] and others for the homogeneous 
materials. 

Consider the relation between D and E directly expressed 
from Faraday's and Ampere's Laws: 
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where ^o is the permeability of free space. Expansion of the 
matrix form of M\\ is 
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Using the central difference operators, Eq. ( I12> can be dis- 
cretized as 
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and 5y and 5^ 



are defined in the same way as in [28]. In 
addition to the wave equation, the constitutive relation of wire 
medium (|8|l must also be considered and can be written in the 
matrix form: 
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For stability analysis in accordance to [31], we substitute the 
following solution into the discrete equations ( I13l l and ( I14> . 

pm _ p^n j(m^A^q^+myAyqy+7n^A,q,) q^x 

where F is a complex amplitude, Z is the complex variable 
which gives the growth of the error in a time interation and 
q = (QxjQyTQz) is the numerical wave vector of the discrete 
mode. After simple calculations we obtain 
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respectively. The determinant of the system of equations ( I17> 
and (^} provides us with the stability polynomial; 
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For simplicity, the terms that will lead to the stability condition 
in free space (i.e. the conventional stability condition) are 
omitted from this stability polynomial. 

In order to avoid numerical root searching [30] for obtaining 
the stability conditions, the above stability polynomial can be 
transformed into the r-plane using the bilinear transformation 

r- -I- 1 

(19) 
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The stability polynomial in the r-plane becomes as follows; 
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Building up the Routh table for the above polynomial as in 
[31], we obtain the following stability conditions; 
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In order to fulfil these conditions, it is enough to fulfil the 
conventional Courant stability condition [6]; 

\a^x,y,z " / 

Therefore the conventional Courant stability condition [6] is 
preserved for modelling of wire medium and no additional 
conditions are required. 

Note that in the above analysis, the central average operator 
/Zj was used for discretization of the Eq. (|6}. If we choose the 
central average operator /i2t defined as 

h^2t^ Irrix^my^mz \ Irn^^rny^mz ' Irrix^my^mz j / ' \-^^J 

then the stability condition ( 12 1> remains the same, but ( 12 2t 
becomes 
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which indicates that even less restrictive stability condition 
than that for the conventional FDTD method can be reached. 
However, if no central average operator is used when discretis- 
ing the Eq. (|6j, then (I21t also remains the same, but (I22t will 
change to 

,2, 
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sm 
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Therefore, if the plasma frequency used in simulations is 
too high and the time step is not properly chosen, then the 
discretised formulations can become unstable. Fig. |3] shows 
the comparison of magnetic field distribution using different 
discretisation schemes; using the central average operator /ij 
and without using central average operator It is clearly shown 
that after 500 time steps, the instability errors start appearing 
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Fig. 3. Comparison of the internal magnetic field distributions (unit: A/m) 
in the plane y = for a wire medium slab excited by a point magnetic 
source (d = 0.5A, w = X, h = O.IA and ko/k = 4) calculated using 
different discretisation schemes: dashed line - the central average operator 
/x^, solid line - without using central average operator The field is plotted 
at the time step n = 570 Af, where At is chosen as the Courant limit i.e. 
At = A^/V2c, A^ = A/200. 



Fig. 4. Reflection en'or (in dB) from the MPML-wire medium and PML- 
free space interfaces calculated at the observation plane (2 cells away from the 
PMLs) for a wire medium slab excited by a point magnetic source (d = 0.5A, 
h = O.IA and ko/k = 4) plotted as functions of the time step. The wire 
medium along {/-direction is long enough to ensure the wave reflected back 
from the far boundary does not reach the reference plane during calculations. 



from inside the wire medium slab using the latter scheme. 

B. Numerical Dispersion 

The numerical dispersion relation for the wire medium 
can be found by evaluating the stability polynomial Sw(Z) 
given by (I18> on the unit circle of the Z-plane (i.e. by 
letting Z — e^^^*), and equating the results to zero. After 
some calculations, the numerical dispersion relation for wire 
medium is obtained: 
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If A^ ^ 0, where (3 = x,y,z,t, then ( I27> reduces to the 
continuous dispersion relation for wire medium [4]: 



{ql 



ill 



Qz 







(28) 



The first and second terms of j28> coiTespond to the transmis- 
sion line modes (TEM waves with respect to the orientation 
of wires) and the extraordinary modes (TM waves), see [4] 
for details. The ordinary modes (TE waves) do not appear in 
(I28> since their contribution was omitted in ( I18> for simplicity 
of calculation. 

IV. Perfectly Matched Layer Formulation 

In 1994, Berenger introduced a nonphysical absorber for 
terminating the outer boundaries of the FDTD computation 
domain that has a wave impedance independent of the angle 
of incidence and frequency. This absorber is called the per- 
fectly matched layer (PML) [32]. The development of PML 
involves a splitting-field approach and the reflection from a 
PML boundary is dependent only on the PML's depth and 
conductivity. In [33], Berenger's original PML is extended 



to absorb electromagnetic waves propagating in anisotropic 
dielectric and magnetic media by introducing the material- 
independent quantities (electric flux density D and magnetic 
flux density B). 

PMLs are usually placed at a distance of A/2 from any 
objects in the simulation domain. In order to reduce the time 
and computer memory requirements for simulations as well 
as to improve the convergence performance, it is required to 
place the PML in the close vicinity of the wire medium. For 
that purpose, we can follow a similar approach as in [33] 
by modifying Berenger's original PML formulations. In the 
modified PML for wire medium, D is introduced into the 
updating equations and the quantities D and H are splitted. 
For example, the updating equation for D^x becomes 
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(29) 



and the updating equations for H remain the same as in 
Berenger's original PML [32]. The matching conditions are, 

c^f = erf , where a = x,y, z, (30) 

where a^ and a^ denote electric conductivity and magnetic 
loss inside the PML, respectively. It should be noted that 
the difference between the above expressions and Berenger's 
original PML formulations is that Eq is not involved in the 
expression of the theoretical reflection factor R{0) for erf 
[33]. Furthermore, the updating relation between D and E in 
Eq. ( I1Q> must be extended into the PML in order to match 
wire medium with the modified PML. 

In order to evaluate the performance of the modified PML 
for the wire medium, a 2-D computation domain similar to that 
in Fig.|2]is chosen except that along y-direction where the wire 
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Fig. 5. (a) Distribution of magnetic field for a finite wire medium slab excited 
by a point magnetic source {d = 0.5A, w = X, h = O.IA and ko/k = 4). 
(b) The same slab, but terminated by a ten-cell MPML at each side along 
j;-direction. (unit: A/m) 



medium slab is directly terminated by a ten-cell modified PML 
(MPML) with a normal theoretical reflection i?(0) = 10"^ 
(see sketch in Fig. 0}. From the other sides the Berenger's 
original PML is used to truncate the free space. The source is 
located close to the edge of the wire medium slab terminated 
by MPML and enough far away to the other side of the slab 
in order to ensure that the wave reflected back from that 
side does not reach the reference plane during calculations. 
The observation plane is 2 cells away from the MPML. The 
reference plane is located at the same distance from the source 
as the observation plane, but from the other side (see sketch in 
Fig. 1^. The magnetic fields at the observation and reference 
planes are recorded as ilpML and Hid, respectively. The 
reflection error is defined as 

Reflection Error (dB) = 20 x logip f ^^^^^~^" 

(31) 
where |i/max| is the maximum value of the magnetic field at 



the reference plane. The reflection error is then calculated for 
Hz and shown in Fig. |4] For comparison, the reflection error 
at the PML-free space interface is also shown. It is found that 
the reflections from the PML are less than -70 dB thus the 
wire medium is "perfectly" matched. 

Fig. |5] shows the comparison of magnetic field distributions 
for two cases; using Berenger's original PML at A/2 distance 
from the slab and using the modified PML to truncate the 
wire medium slab directly. In comparison with Fig. |5ja), the 
simulation domain size for Fig. (Sjb) is reduced by 50% and 
the convergence is greatly improved since the diffractions from 
the corners and edges are avoided in simulations. For the first 
case, the reflection error falls below -30 dB after 1000 periods 
(400,000 time steps), while for the latter one, the convergence 
is reached after 100 periods. 

V. Modelling of The Sub-wavelength Imaging 

In order to validate the proposed spatially dispersive FDTD 
formulations, we have chosen flat sub-wavelength lenses 
formed by wire media [29], [34]. Such lenses provide unique 
opportunity to transfer images with resolution below classical 
diffraction limit. In the present paper we have considered a 
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Fig. 7. Comparison of magnetic field distribution at the source and image 
planes for a wire medium slab (d = 0.5A, h = 0.05A and kg/k = 4) directly 
terminated by a ten-cell MPML at each side along ^/-direction, (unit: A/m) 



2-D case: the structure is infinite in z-direction and electric 
field is in x-y plane (TM polarization with respect to the 
orientation of wires). The following parameters are used in 
simulations: the operating frequency is 3.0 GHz (wavelength 
A = 0.1 m in free space) and the plasma frequency of the wire 
medium is /o = 12.0 GHz (ko/k = 4); the FDTD cefl size is 
Ax = A/2GG with the time step At = 8.33 x 10"" s according 
to the stability criteria ( I23l l: a ten-cell Berenger's original 
PML is used to truncate the free space and the modified PML 
(MPML) is used to match the wire medium slab; the thickness 
of the wire media slab is chosen as d — A/2. Three equally 
spaced (A/20) magnetic point sources with phase differences 
of 180° between neighbouring sources are located at a distance 
of A/20 from the left interface of the wire medium slab. 

Fig. |6l shows the distribution of the x- and y- components 
of electric field as well as z- component of magnetic field in 
the simulation domain. One can see from Fig. |6ja) that x- 
component of electrical field penetrates into the wire medium 
in the form of extraordinary modes [4] which are evanescent 
and decay with distance. That is why this component vanishes 
in the center of the slab. The extraordinary modes are coupled 
with transmission line modes of wire medium which are 
clearly seen in Figs. |6jb) and |6jc). In accordance to the 
canalisation principle [29], the transmission line modes deliver 
image from the front interface to the back one. Fig. shows 
the magnetic field distribution at the source and image planes 
(located on the different sides of wire medium slab as in Fig. 
13 . It is worth noting that the image is in phase or out-of- 
phase with the source if the thickness of wire medium slab is 
even or odd integer numbers of A/2, respectively [29]. That is 
why in the case under consideration the image appears in out 
of phase. It can be seen that in Fig. the distance between 
two maxima is approximately A/IG which verifies the sub- 
wavelength imaging capability of the wire medium lenses. The 
performed simulations confirm that the spatial dispersion in 
wire medium is accurately taken into account in the presented 
spatially dispersive FDTD model. 



Y. ZHAO et al: MODELLING OF WAVE PROPAGATION IN WIRE MEDIA USING FDTD 



a) 



1.8 



?s 



0.8 



0.6 



0.4 



0.2 



0.5 



xlX 




Fig. 6. The distributions of a) Ex (unit: V/m), b) Ey (unit: V/m) and c) H^ (unit: A/m) for a wire medium slab excited by three equally spaced magnetic 
sources with the phase differences equal to tt (d = 0.5A, h = 0.05A and k^/k = 4) directly terminated by a ten-cell MPML at each side along ^/-direction. 



VI. Conclusion 

The spatially dispersive FDTD formulations have been 
developed for modelling of wave propagation in wire medium 
as effective dielectrics. The auxiliary differential equation 
method is used in order to take into account both the spatial 
and frequency dispersion effects. The stability analysis shows 
that the conventional Courant stability limit is preserved if 
the standard central difference approximations and central 
average operator are used to discretise differential equations. 
Through the use of modified PML, the wire medium can 
be "perfectly" matched to the absorbing boundaries thus the 
convergence performance in simulations is greatly improved 
since the diffractions from corners and edges of the finite sized 
wire medium slab are avoided. The flat sub-wavelength lenses 
formed by wire medium are chosen for the validation of devel- 
oped spatially dispersive FDTD formulations. Numerical simu- 
lation results verifies the sub-wavelength imaging capability of 
wire media. The proposed spatially dispersive FDTD method 
will demonstrate its distinct value when complex sources are 
used in the simulation for exploitation of practical applications 
of wire medium in antenna and microwave engineering. 
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